//！ 这个基础数学计算是不是应该把返回的值改为数值
//！ 因为这样的话，容易并行计算
//! 处理错误的话不容易并行计算
//! 或者对于并行计算库比如futhark,再重写代码

const std = @import("std");
const math = std.math;
const dprint = std.debug.print;

// pub fn LB_Math(comptime T: type) type {
//     return struct{

//     };
// }

pub const Math_Compute_Abandon = error{
    Math_Compute_Abandon,
};

test defaultTolerance {
    const tol_f64 = defaultTolerance(f64);
    const tol_f32 = defaultTolerance(f32);
    std.debug.print("default tolerace :{} {}\n", .{ tol_f64, tol_f32 });
}

/// 浮点类型相关的默认容差
fn defaultTolerance(comptime T: type) T {
    return switch (@typeInfo(T)) {
        .float, .comptime_float => math.floatEps(T) * 100, // 100倍机器精度
        else => @compileError("Only floating-point types supported"),
    };
}

/// 根据绝对误差判断是否大致相等
pub fn approxEqAbs(
    comptime T: type,
    x: T,
    y: T,
    tolerance: ?T, // 可选参数，为null时使用类型默认值
) bool {
    if (x == y) {
        return true;
    }
    //     // 处理NaN和无穷大
    if (math.isNan(x) or math.isNan(y)) {
        return false;
    }
    const tol = tolerance orelse defaultTolerance(T);
    return @abs(x - y) <= tol;
}

test safe_sqrt {
    var x: f32 = -1;
    x = x - 1;
    try std.testing.expect(safe_sqrt(x) >= 0);
    const n: u32 = 4;
    std.debug.print("factorial :{}\n", .{factorial(n)});
    // math.divExact(comptime T: type, numerator: T, denominator: T)
}

fn Sqrt(comptime T: type) type {
    return switch (@typeInfo(T)) {
        .int => |int| @Type(.{ .int = .{ .signedness = .unsigned, .bits = (int.bits + 1) / 2 } }),
        else => T,
    };
}

pub fn safe_sqrt(x: anytype) Sqrt(@TypeOf(x)) {
    if (x < 0) {
        return math.sqrt(0);
    }
    return math.sqrt(x);
}

pub fn factorial(n: u64) u128 {
    var re: u128 = 1;
    var i: u128 = 2;
    while (i <= n) : (i += 1) {
        re *= i;
    }
    return re;
}

pub fn pailie_map_i_to_reverse_order(n: u64, arr: []u64) void {
    const num = factorial(n);
    // std.debug.assert(arr.len == n * num);

    for (0..num) |i| {
        var chu = i;
        for (0..n) |j| {
            const yu = chu % (j + 1);
            chu = chu / (j + 1);
            arr[i * n + n - 1 - j] = yu;
        }
    }
}

// test map_reverse_order_to_pailie {}
pub fn map_reverse_order_to_pailie(arr: []u64, n: u64, re: []u64) void {
    const num = factorial(n);
    for (0..num) |i| {
        for (0..n) |j| {
            re[i * n + j] = j;
        }
    }
    // var ix = 0;
    // var iy = 0;
    // var temp = 0;
    for (0..num) |i| {
        for (0..n) |j| {
            const ix = i * n + j;
            const iy = i * n + j + arr[i * n + j];
            var l = iy;
            while (l > ix) : (l -= 1) {
                const temp = re[l - 1];
                re[l - 1] = re[l];
                re[l] = temp;
            }
        }
    }
}
fn increase_a_step_for_combination(m: u64, n: u64, temp_array: []u64) Math_Compute_Abandon!void {
    var mu_arr: []u64 = temp_array;
    var temp_i: u64 = 0;

    for (0..m) |i| {
        temp_i = mu_arr[i] + 1;
        if (temp_i < n) {
            if ((i + 1) < m) {
                if (temp_i >= mu_arr[i + 1]) {
                    continue;
                }
            }
            // var temp_mu = mu_arr[m..];
            // std.mem.copyForwards(comptime T: type, dest: []T, source: []const T)
            for (0..m) |j| {
                mu_arr[j + m] = mu_arr[j];
            }
            mu_arr[i + m] = temp_i;
            // mu_arr = temp_mu;
            for (0..i) |j| {
                mu_arr[j + m] = j;
            }
            return;
        } else {
            return error.Math_Compute_Abandon;
        }
    }
    return error.Math_Compute_Abandon;
}

pub fn compute_combination(m: u64, n: u64, re: []u64) Math_Compute_Abandon!void {
    if (m > n or m <= 0 or n <= 0) {
        return error.Math_Compute_Abandon;
    }
    var re_len: u128 = 1;
    for (0..m) |i| {
        re_len *= (n - i);
    }
    re_len = re_len / factorial(m);
    for (0..m) |i| {
        re[i] = i;
    }
    if (re.len < re_len * m) {
        return error.Math_Compute_Abandon;
    }
    for (0..(re_len - 1)) |i| {
        try increase_a_step_for_combination(m, n, re[i * m .. (i + 2) * m]);
    }
}

// test lb_norm {
//     var p: [3]f64 = .{ 1, 2, 3 };
//     p[0] = 2;
//     const re = lb_norm(f64, &p);
//     dprint("lb norm :{}\n", .{re});
// }

test lb_norm {
    var p: [3]f64 = .{ 1, 2, 3 };
    p[0] = 2;
    const re = lb_norm(&p);
    dprint("lb norm :{}\n", .{re});
}
pub fn lb_norm(p: anytype) Sqrt(@TypeOf(p[0])) {
    var re: @TypeOf(p[0]) = 0;
    for (p) |v| {
        re += v * v;
    }
    return safe_sqrt(re);
}

test distance_of_two_points {
    const a: [4]f64 = .{ 1, 2, 3, 4 };
    const b: [4]f64 = .{ 4, 3, 2, 1 };
    dprint("distance_of_two_points:{}\n", .{distance_of_two_points(&a, &b)});
}
pub fn distance_of_two_points(a: anytype, b: []const (@TypeOf(a[0]))) Sqrt(@TypeOf(a[0])) {
    var re: @TypeOf(a[0]) = 0;
    for (a, b) |va, vb| {
        re += (va - vb) * (va - vb);
    }
    return safe_sqrt(re);
}

// test out_product {
//     var a = [3]f64{ 1.0, 2.0, 3.0 };
//     var b = [3]f64{ 0.0, 0.5, -2.0 };
//     a[0] = a[1] * b[2];
//     b[1] = a[2] - b[0];

//     const out = out_product(f64, &a, &b);
//     std.debug.print("out :{any}\n", .{out});
// }

test out_product {
    var a = [3]f64{ 1.0, 2.0, 3.0 };
    var b = [3]f64{ 0.0, 0.5, -2.0 };
    a[0] = a[1] * b[2];
    b[1] = a[2] - b[0];

    const out = out_product(&a, &b);
    std.debug.print("out :{any}\n", .{out});
}

test "comptiem typeinfo" {
    const a: [10]f64 = [_]f64{0} ** 10;
    const b = &a;
    var len: usize = 2;
    len *= 2;
    len *= 2;
    const c = a[0..len];

    std.debug.print("comptime typeinfo tyep {}\n", .{@typeInfo(@TypeOf(b)).pointer.size});
    std.debug.print("comptime typeinfo1 tyep {}\n", .{@typeInfo(@TypeOf(c)).pointer.size});
}

// 调用时最好直接两个变量直接切片[0..3]
pub fn out_product(a: anytype, b: *const [3](@TypeOf(a[0]))) [3]@TypeOf(a[0]) {
    // std.debug.assert(a.len == 3 and b.len == 3 and re.len == 3);
    comptime {
        const P = @TypeOf(a);
        switch (@typeInfo(P)) {
            .pointer => |p| {
                switch (p.size) {
                    .slice => if (a.len != 3) @compileError("a.len is not 3\n"),
                    .one => if (p.child != [3]@TypeOf(a[0])) @compileError("type erro"),
                    else => @compileError("type erro\n"),
                }
            },
            .array => if (a.len != 3) @compileError("a.len is not 3\n"),
            else => @compileError("not pointer\n"),
        }
    }
    var re: [3]@TypeOf(a[0]) = .{ 0, 0, 0 };
    re[0] = a[1] * b[2] - a[2] * b[1];
    re[1] = -a[0] * b[2] + a[2] * b[0];
    re[2] = a[0] * b[1] - a[1] * b[0];
    return re;
}

pub fn inner_product(a: anytype, b: []const @TypeOf(a[0])) @TypeOf(a[0]) {
    comptime {}
    var re: @TypeOf(a[0]) = 0;
    for (a, b) |va, vb| {
        re += va * vb;
    }
    return re;
}

test normalize {
    var arr1: [3]f32 = .{ -1.3, 2.1, -6.2 };
    try normalize(&arr1);
    dprint("{any} {}\n", .{ arr1, inner_product(arr1, &arr1) });
}
pub fn normalize(p: anytype) Math_Compute_Abandon!void {
    comptime {
        const T = @TypeOf(p[0]); // 推导元素类型 T
        const P = @TypeOf(p); // p 的完整类型
        const info = @typeInfo(P);

        // 检查是否为可变切片 []T
        const is_mutable_slice =
            info == .pointer and
            info.pointer.size == .slice and
            !info.pointer.is_const and // 关键：使用 is_const 字段
            info.pointer.child == T;

        // 检查是否为指向数组的可变指针 *[N]T
        const is_ptr_to_array =
            info == .pointer and
            info.pointer.size == .one and
            !info.pointer.is_const and // 关键：指针本身不能是 const
            @typeInfo(info.pointer.child) == .array;

        // 如果两种都不符合，编译失败
        if (!is_mutable_slice and !is_ptr_to_array) {
            @compileError(
                \\normalize() expects a mutable slice `[]{s}` or a pointer to an array `*[{d}]{s}`,
                \\but got `{s}`.
                \\
                \\Note: `[]const {s}` and `*const [{d}]{s}` are not allowed because the function modifies the data.
            ++ @typeName(T) ++ @typeName(T) ++ @typeName(P) ++
                @typeName(T) ++ @typeName(T));
        }
    }
    const norm = lb_norm(p);
    if (approxEqAbs(@TypeOf(p[0]), norm, 0, null)) {
        return error.Math_Compute_Abandon;
    }
    for (p) |*pv| {
        pv.* /= norm;
    }
}
// tent映射，均匀分布

pub fn tent_map(x: anytype) @TypeOf(x) {
    return if (x < 0.5) 2.0 * x else 2.0 * (1 - x);
    // return x<0.5?(2.0*x):(2.0*(1-x));
}
// sine map
pub fn sin_map(x: anytype) @TypeOf(x) {
    // math.sin
    return math.sin(math.pi * x);
}

fn quaternions_plus(comptime T: type, q1: [4]T, q2: [4]T) [4]T {
    return .{
        q1.v[0] + q2.v[0],
        q1.v[1] + q2.v[1],
        q1.v[2] + q2.v[2],
        q1.v[3] + q2.v[3],
    };
}

fn quaternions_mult(comptime T: type, q1: [4]T, q2: [4]T) [4]T {
    var re: [4]T = [_]T{0} ** 4;
    re[0] = q1[0] * q2[0] - q1[1] * q2[1] - q1[2] * q2[2] - q1[3] * q2[3];

    re[1] = q1[0] * q2[1] + q1[1] * q2[0] + q1[2] * q2[3] - q1[3] * q2[2];

    re[2] = q1[0] * q2[2] + q1[2] * q2[0] + q1[3] * q2[1] - q1[1] * q2[3];

    re[3] = q1[0] * q2[3] + q1[3] * q2[0] + q1[1] * q2[2] - q1[2] * q2[1];

    return re;
}

fn quaternions_conjugate(comptime T: type, q: [4]T) [4]T {
    return .{
        q[0],
        -q[1],
        -q[2],
        -q[3],
    };
}
fn q_rota(comptime T: type, q1: [4]T, x: [3]T) [4]T {
    const re = [4]T{ 0, x[0], x[1], x[2] };

    const q_ = quaternions_conjugate(T, q1);

    return quaternions_mult(T, q1, quaternions_mult(T, re, q_));
}

pub fn quaternions_inverse(comptime T: type, q: [4]T) [4]T {
    // const norm = lb_norm(T, q);
    const norm = inner_product(q, &q);
    return .{ q[0] / norm, -q[1] / norm, -q[2] / norm, -q[3] / norm };
}

pub fn get_rotation_from_matrix(comptime T: type, matrix: *const [3][3]T) [4]T {
    var re = [_]T{0} ** 4;
    // double* re=(double*)malloc(sizeof(double)*4);
    // memset(re,0,sizeof(double)*4);
    re[0] = math.acos((matrix[0][0] + matrix[1][1] + matrix[2][2] - 1) / 2.0);

    // re[0]=acos( (matrix[0]+matrix[4]+matrix[8]  -1)/2  );

    // if(fabs(re[0])<1e-9)
    if (approxEqAbs(T, re[0], 0, null)) {
        re[0] = 0;
        re[1] = 1.0;
        return re;
    }
    re[1] = (matrix[2][1] - matrix[1][2]) / (2.0 * math.sin(re[0]));
    re[2] = (matrix[0][2] - matrix[2][0]) / (2.0 * math.sin(re[0]));
    re[3] = (matrix[1][0] - matrix[0][1]) / (2.0 * math.sin(re[0]));
    return re;
}

// double*  get_rotation_from_quaternions( Quaternions q  )
pub fn get_rotation_from_quaternions(comptime T: type, q: [4]T) [4]T {

    // double* re=(double*)malloc(sizeof(double)*4);
    var re = [_]T{0} ** 4;

    re[0] = 2 * math.acos(q[0]);
    re[1] = q[1];
    re[2] = q[2];
    re[3] = q[3];
    normalize(re[1..]);
    return re;
}
// static inline double* get_matrix_from_rotation(double theta,double* axis)

pub fn get_matrix_from_rotation(comptime T: type, theta: T, axis: [3]T) [3][3]T {
    // double *re=(double*)malloc(sizeof(double)*9);
    var re: [3][3]T = [3][3]T{
        [3]T{ 1, 0, 0 },
        [3]T{ 0, 1, 0 },
        [3]T{ 0, 0, 1 },
    };
    normalize(&axis) catch return re;

    const sint = math.sin(theta);
    const cost = math.cos(theta);
    // double sint=sin(theta);
    // double cost=cos(theta);
    re[0][0] = cost + (1 - cost) * axis[0] * axis[0];
    re[0][1] = (1 - cost) * axis[0] * axis[1] - axis[2] * sint;
    re[0][2] = (1 - cost) * axis[0] * axis[2] + axis[1] * sint;

    re[1][0] = sint * axis[2] + (1 - cost) * axis[0] * axis[1];
    re[1][1] = cost + (1 - cost) * axis[1] * axis[1];
    re[1][2] = (1 - cost) * axis[1] * axis[2] - axis[0] * sint;

    re[2][0] = (1 - cost) * axis[0] * axis[2] - axis[1] * sint;
    re[2][1] = (1 - cost) * axis[2] * axis[1] + axis[0] * sint;
    re[2][2] = (1 - cost) * axis[2] * axis[2] + cost;

    return re;
}

// double*  get_matrix_from_quaternions( Quaternions q  )
fn get_matrix_from_quaternions(comptime T: type, q: [4]T) [3][3]T {
    const theta = 2 * math.acos(q[0]);
    var ax: [3]T = .{ q[1], q[2], q[3] };
    normalize(&ax);
    // double ax[3]={q.v[1],q.v[2],q.v[3]};
    // normalize(ax,3);
    return get_matrix_from_rotation(T, theta, ax);
}

// static inline Quaternions quaternions_power(Quaternions q,double t)
pub fn quaternions_power(comptime T: type, q: [4]T, t: T) [4]T {
    // Quaternions re;
    // quaternions_init(&re);
    var re: [4]T = [_]T{0} ** 4;
    const norm = lb_norm(&q);
    const theta = math.acos(q[0] / norm);
    // double norm= quaternions_norm(q);
    // double theta=acos(q.v[0]/norm);
    const a = math.pow(norm, t);
    // double a=pow(norm,t);
    re[0] = a * math.cos(theta * t);
    re[1] = a * math.sin(theta * t) * (q[1] / (norm * math.sin(theta)));
    re[2] = a * math.sin(theta * t) * (q[2] / (norm * math.sin(theta)));
    re[3] = a * math.sin(theta * t) * (q[3] / (norm * math.sin(theta)));

    return re;
}

pub fn q_rotas_interpolation(comptime T: type, q1: [4]T, q2: [4]T, t: T) [4]T {
    return quaternions_mult(T, quaternions_power(T, q1, 1 - t), quaternions_power(T, q2, t));
}
// 四元数测试代码
// test "quaternions basic operations" {
//     const T = f64;

//     // 测试四元数乘法
//     const q1 = [4]T{ 1.0, 2.0, 3.0, 4.0 };
//     const q2 = [4]T{ 5.0, 6.0, 7.0, 8.0 };
//     const product = quaternions_mult(T, q1, q2);

//     // 验证四元数乘法结果 (手工计算验证)
//     try std.testing.expectApproxEqAbs(@as(T, -60.0), product[0], 1e-10);
//     try std.testing.expectApproxEqAbs(@as(T, 12.0), product[1], 1e-10);
//     try std.testing.expectApproxEqAbs(@as(T, 30.0), product[2], 1e-10);
//     try std.testing.expectApproxEqAbs(@as(T, 24.0), product[3], 1e-10);

//     // 测试共轭
//     const conjugate = quaternions_conjugate(T, q1);
//     try std.testing.expectApproxEqAbs(@as(T, 1.0), conjugate[0], 1e-10);
//     try std.testing.expectApproxEqAbs(@as(T, -2.0), conjugate[1], 1e-10);
//     try std.testing.expectApproxEqAbs(@as(T, -3.0), conjugate[2], 1e-10);
//     try std.testing.expectApproxEqAbs(@as(T, -4.0), conjugate[3], 1e-10);

//     // 测试单位四元数
//     const identity = [4]T{ 1.0, 0.0, 0.0, 0.0 };
//     const product_with_identity = quaternions_mult(T, q1, identity);
//     try std.testing.expectApproxEqAbs(q1[0], product_with_identity[0], 1e-10);
//     try std.testing.expectApproxEqAbs(q1[1], product_with_identity[1], 1e-10);
//     try std.testing.expectApproxEqAbs(q1[2], product_with_identity[2], 1e-10);
//     try std.testing.expectApproxEqAbs(q1[3], product_with_identity[3], 1e-10);
// }

// test "quaternions rotation" {
//     const T = f64;

//     // 创建绕Z轴旋转90度的四元数
//     const angle = math.pi / 2.0;
//     const axis = [3]T{ 0.0, 0.0, 1.0 };
//     const q = [4]T{ math.cos(angle / 2.0), axis[0] * math.sin(angle / 2.0), axis[1] * math.sin(angle / 2.0), axis[2] * math.sin(angle / 2.0) };

//     // 测试点(1, 0, 0)旋转90度后应该变为(0, 1, 0)
//     const point = [3]T{ 1.0, 0.0, 0.0 };
//     const rotated = q_rota(T, q, point);

//     // 旋转后的点应该是纯四元数，实部为0
//     try std.testing.expectApproxEqAbs(@as(T, 0.0), rotated[0], 1e-10);
//     try std.testing.expectApproxEqAbs(@as(T, 0.0), rotated[1], 1e-10);
//     try std.testing.expectApproxEqAbs(@as(T, 1.0), rotated[2], 1e-10);
//     try std.testing.expectApproxEqAbs(@as(T, 0.0), rotated[3], 1e-10);
// }

test "quaternions inverse" {
    const T = f64;

    // 创建一个四元数
    const q = [4]T{ 1.0, 2.0, 3.0, 4.0 };
    const inv_q = quaternions_inverse(T, q);

    // 四元数与其逆的乘积应该是单位四元数
    const product = quaternions_mult(T, q, inv_q);

    // 注意：对于非单位四元数，乘积的实部不为1，但虚部应该为0
    try std.testing.expectApproxEqAbs(@as(T, 1.0), product[0], 1e-10);
    try std.testing.expectApproxEqAbs(@as(T, 0.0), product[1], 1e-10);
    try std.testing.expectApproxEqAbs(@as(T, 0.0), product[2], 1e-10);
    try std.testing.expectApproxEqAbs(@as(T, 0.0), product[3], 1e-10);
}

// test "rotation matrix conversions" {
//     const T = f64;

//     // 创建绕Z轴旋转90度的旋转矩阵
//     const angle = math.pi / 2.0;
//     const axis = [3]T{ 0.0, 0.0, 1.0 };
//     const matrix = get_matrix_from_rotation(T, angle, axis);

//     // 从矩阵提取旋转四元数
//     const q_from_matrix = get_rotation_from_matrix(T, &matrix);

//     // 从四元数提取旋转轴和角度
//     const rotation = get_rotation_from_quaternions(T, q_from_matrix);

//     // 验证旋转角度
//     try std.testing.expectApproxEqAbs(angle, rotation[0], 1e-10);

//     // 验证旋转轴 (应该是Z轴)
//     const norm = math.sqrt(rotation[1] * rotation[1] + rotation[2] * rotation[2] + rotation[3] * rotation[3]);
//     try std.testing.expectApproxEqAbs(axis[0], rotation[1] / norm, 1e-10);
//     try std.testing.expectApproxEqAbs(axis[1], rotation[2] / norm, 1e-10);
//     try std.testing.expectApproxEqAbs(axis[2], rotation[3] / norm, 1e-10);

//     // 从四元数转换回矩阵
//     const matrix_from_q = get_matrix_from_quaternions(T, q_from_matrix);

//     // 验证矩阵是否相同
//     for (0..3) |i| {
//         for (0..3) |j| {
//             try std.testing.expectApproxEqAbs(matrix[i][j], matrix_from_q[i][j], 1e-10);
//         }
//     }
// }

// test "quaternions power and interpolation" {
//     const T = f64;

//     // 创建两个单位四元数
//     const q1 = [4]T{ math.cos(math.pi / 4.0), 0.0, 0.0, math.sin(math.pi / 4.0) }; // 绕Z轴旋转90度
//     const q2 = [4]T{ math.cos(math.pi / 2.0), 0.0, 0.0, math.sin(math.pi / 2.0) }; // 绕Z轴旋转180度

//     // 测试四元数幂运算
//     const q_half = quaternions_power(T, q2, 0.5);

//     // q2的平方根应该等于q1 (因为180度的平方根是90度)
//     try std.testing.expectApproxEqAbs(q1[0], q_half[0], 1e-10);
//     try std.testing.expectApproxEqAbs(q1[1], q_half[1], 1e-10);
//     try std.testing.expectApproxEqAbs(q1[2], q_half[2], 1e-10);
//     try std.testing.expectApproxEqAbs(q1[3], q_half[3], 1e-10);

//     // 测试四元数插值
//     const t = 0.5;
//     const q_interpolated = q_rotas_interpolation(T, q1, q2, t);

//     // 在中间点插值应该得到135度的旋转
//     const expected_angle = 3.0 * math.pi / 4.0;
//     const expected_q = [4]T{ math.cos(expected_angle / 2.0), 0.0, 0.0, math.sin(expected_angle / 2.0) };

//     try std.testing.expectApproxEqAbs(expected_q[0], q_interpolated[0], 1e-10);
//     try std.testing.expectApproxEqAbs(expected_q[1], q_interpolated[1], 1e-10);
//     try std.testing.expectApproxEqAbs(expected_q[2], q_interpolated[2], 1e-10);
//     try std.testing.expectApproxEqAbs(expected_q[3], q_interpolated[3], 1e-10);
// }

// test "quaternions normalization" {
//     const T = f64;

//     // 创建非单位四元数
//     var q = [4]T{ 2.0, 4.0, 6.0, 8.0 };

//     // 归一化
//     const norm = math.sqrt(q[0] * q[0] + q[1] * q[1] + q[2] * q[2] + q[3] * q[3]);
//     q[0] /= norm;
//     q[1] /= norm;
//     q[2] /= norm;
//     q[3] /= norm;

//     // 验证归一化后的四元数是单位四元数
//     const product = quaternions_mult(T, q, quaternions_conjugate(T, q));

//     try std.testing.expectApproxEqAbs(@as(T, 1.0), product[0], 1e-10);
//     try std.testing.expectApproxEqAbs(@as(T, 0.0), product[1], 1e-10);
//     try std.testing.expectApproxEqAbs(@as(T, 0.0), product[2], 1e-10);
//     try std.testing.expectApproxEqAbs(@as(T, 0.0), product[3], 1e-10);
// }

// // 测试四元数加法
// test "quaternions addition" {
//     const T = f64;

//     const q1 = [4]T{ 1.0, 2.0, 3.0, 4.0 };
//     const q2 = [4]T{ 5.0, 6.0, 7.0, 8.0 };
//     const sum = quaternions_plus(T, q1, q2);

//     try std.testing.expectApproxEqAbs(@as(T, 6.0), sum[0], 1e-10);
//     try std.testing.expectApproxEqAbs(@as(T, 8.0), sum[1], 1e-10);
//     try std.testing.expectApproxEqAbs(@as(T, 10.0), sum[2], 1e-10);
//     try std.testing.expectApproxEqAbs(@as(T, 12.0), sum[3], 1e-10);
// }

// // 测试四元数旋转向量的应用
// test "quaternion vector rotation" {
//     const T = f64;

//     // 创建绕Y轴旋转180度的四元数
//     const angle = math.pi;
//     const axis = [3]T{ 0.0, 1.0, 0.0 };
//     const q = [4]T{ math.cos(angle / 2.0), axis[0] * math.sin(angle / 2.0), axis[1] * math.sin(angle / 2.0), axis[2] * math.sin(angle / 2.0) };

//     // 测试点(1, 0, 0)旋转180度后应该变为(-1, 0, 0)
//     const point = [3]T{ 1.0, 0.0, 0.0 };
//     const rotated = q_rota(T, q, point);

//     // 旋转后的点应该是纯四元数，实部为0
//     try std.testing.expectApproxEqAbs(@as(T, 0.0), rotated[0], 1e-10);
//     try std.testing.expectApproxEqAbs(@as(T, -1.0), rotated[1], 1e-10);
//     try std.testing.expectApproxEqAbs(@as(T, 0.0), rotated[2], 1e-10);
//     try std.testing.expectApproxEqAbs(@as(T, 0.0), rotated[3], 1e-10);
// }
